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A mathematical model for phase separation: 
a generalized Cahn-Hilliard equation. 

A. Berti 3 and I. Bochicchio* 

In this paper we present a mathematical model to describe the phenomenon of phase separation, which is modelled as 
space regions where an order parameter changes smoothly. The model proposed, including thermal and mixing effects, is 
deduced for an incompressible fluid, so the resulting differential system couples a generalized Cahn-Hilliard equation with 
<^> the Navier-Stokes equation. Its consistency with the second law of thermodynamics in the classical Clausius-Duhem form 
is finally proved. 
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1. Introduction 



The mechanism by which a mixture of two or more components separate into distinct regions (or phases) with different chemical 
compositions and physical properties is usually named spinodal decomposition or phase separation. This mechanism differs from 
classical nucleation in that phase separation is much more subtle, and occurs uniformly throughout the material, not just at 
' discrete nucleation sites. 

in 

' Typically the phenomenon of spinodal decomposition occurs when a mixture of two different species, say A and 8, forming a 
single homogeneous phase at a temperature 9 m greater than the critical temperature do, is rapidly cooled to a temperature where 
the homogeneous state is unstable. The resulting inherent instability leads to composition fluctuations, and thus to instantaneous 
phase separation. 

The most common experimental examples of spinodal decomposition occur in metallic alloys [6, 29] and glassy mixtures 
[2, 30]. For example an Al-rich Al-Zn alloy, when quenched rapidly from above 400 °C and then annealed at temperatures in 
the neighborhood of 100 °C, is known to decompose into Al- and Zn-rich regions via the spinodal mechanism [24]. 

The basic theory of spinodal decomposition has been developed, primarily from a metallurgic point of view by Hillert [21], 
Cahn [8, 9], Hilliard [22] and Cook [10]. Subsequently Cahn developed a more general linearized theory of spinodal instability 
pointing out the essential role played by nonlinear effects in determining the nature of the instability and then in limiting its 
growth [24]. 

The phase separation is often described in the framework of phase-field modelling, in that the interface between the two 
pure phases is not sharp but is regarded as a region of finite width having a gradual variation of different physical quantities. 
In addition, to distinguish one phase from the other, it is necessary to select a quantity which differs in the two phases. Since 
Landau, such a quantity is called order parameter and it assumes distinct values in the bulk phases away from the interfacial 
regions over which it varies smoothly. 
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Interpreting the order parameter as the concentration of one of the two metallic components of the binary alloy, Cahn and 
Hilliard [7, 8] introduce the so-called Cahn-Hilliard equation which describes the evolution of the concentration field in a binary 
alloy. 

In the present paper, we present a generalized mathematical model capable of describing a phase separation into the Cahn- 
Hillard theory. Precisely, we consider a mixture of two incompressible fluids with comparable densities but different viscosity, 
and we assume that our system can be described by a single scalar order parameter c, which we can visualize as the difference 
of local mass fraction (concentration) of the two components of a binary solution. In addition, we suppose that the density of 
mixture does not depend on the composition of the mixture {i.e. p(c) = po) such that the general mass balance equation of 
the mixture degenerates to the solenoidal condition. Following [25], in Sect. 3, we postulate that c obeys a diffusion equation. 

The aim of our paper is to propose a model accounting for the fluid motion. In particular, besides the classical coupling 
between the Cahn-Hilliard and Navier-Stokes equations, due to the presence of the material derivative of the order parameter 
in the Cahn-Hilliard equation and of a surface tension source term in the Navier-Stokes equation (see [19, 25]), we suppose 
that the chemical potential may depend on the curl of the velocity of the mixture (see Sect. 4). The effect of the velocity can 
be interpreted as an increase of the temperature which controls the phase separation. 

In Sect. 5, we modify the classical Navier-Stokes equation by adding a reactive stress, accounting for the capillary forces due 
to surface tension, and a skew tensor consistent with the presence of internal structure due to the mixture [11] which guarantees 
the coupling with the Cahn-Hilliard equation. Finally, in Sect. 6, we prove the compatibility with thermodynamics of our model 
by expressing the second law in the classical Clausius-Duhem inequality. 



2. Phase-field modelling 

Let (lei 3 be a fixed bounded domain, which is completely filled by a mixture of two incompressible fluids A and 6, and let 
dQ be its smooth boundary with unit outward normal n. For the sake of simplicity, we suppose that the densities pA, Pb of both 
components as well as the density po of the mixture are constant and we assume 

Pa = Pb = Po = const. (1) 

Let m be the total mass of the mixture, i.e. 

m = p dv 
Jn 

and let m A , m B be the masses of each species in fi, so that m = iva + m B - We denote by Pa, Pb the apparent densities of A 
and B respectively, namely 

m A = pAdv, m B = / pBdv. 
Jn Jn 

As a consequence, 

PO = pA+pB- (2) 

During phase separation each material particle cannot change its phase, but it is only allowed to migrate from a geometrical 
point to another close to it. As a consequence, the total amount of each species in the whole domain must remain equal to the 
given original amount. 

In our model we consider the mixture as a single fluid obeying the laws of conservation of mass and linear momentum of 
continuum mechanics and we associate to each particle of the matter an additional scalar function c, called order parameter, 
which allows us to distinguish one phase (fluid form) from the other one. More precisely, we let c = — 1 in regions filled only by 
the fluid A and c = 1 in regions where only the fluid B appears. 

Following the phase-field approach, we suppose that the two immiscible fluids are not separated by a sharp interface, but we 
assume that there exists a partial mixing between them in thin layers with finite thickness called diffuse interfaces. Accordingly, c 
does not take its values only in { — 1, 1}, but it is allowed to vary smoothly between —1 and 1 in the interfacial regions. Moreover, 
if we suppose rriA = rn B , the condition c = means that the fluid is in a uniform mixed state. Such an approach traces back 
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to van der Waals, Landau and Ginzburg, Cahn and Hilliard ([7, 23, 31]) and later it has been developed in the theory of phase 
transitions (see [1, 4, 5, 14, 19, 27, 28] and the references therein). 

The function c, which we attach to each particle like as a label, may be interpreted as the difference of local mass fraction 
(concentration) of the two components, that is 

c dm = drriA — dm B 

or equivalently 

c _ Pa - Pb 
Po 

In this way, it is apparent from equality (2) that c e [—1, 1]. In particular, c = — 1 (or c = 1) wherever only the component A 
(or 6) occurs. 

Furthermore, the definition of c guarantees that the concentration difference of the two components is conserved in ft as the 
system evolves. Indeed, recalling the definition of pA and pe, we have 

/ poc dv = I (pA — Ps)dv = rriA — ms = constant. 
Jq Jq. 

In the next section, we exhibit a kinetic equation for c able to guarantee the conservation of the total concentration over the 
whole domain. 



Remark 1 Several authors (see [19, 25], for instance) interpret the order parameter c as the local concentration of one of the 
component of the binary fluid. In our notation, it means that c is defined as 



Pa 
Po ' 



As a consequence, c 6 [0, 1]. 



Henceforth, we denote by t the time variable, x, v the position vector and the velocity of the particle at time t in the actual 
configuration, 9 the absolute temperature. Also, V is the gradient operator, the superposed dot is the material derivative and d x 
denotes the partial derivative with respect to the variable %■ In particular, d t is the partial time derivative and 9/ = d Xj . Hence, 
for any function g(x, £), we have 

g = d t g + v-Vg, (3) 

where ■ stands for the scalar product. In addition, we use the symbols V- and A to indicate the divergence and the Laplacian 
respectively. Finally, the inner product of two second order tensors A and B is defined by 

A : B = tr(A T B), 

where trA and A T are the trace and the transpose of a tensor A. 

For reader's convenience, we briefly recall a lemma we will use in the sequel. 



Lemma 1 For any C 2 function g(x, t) the derivatives Vg and Vg are related by the identity 

Wg~ = Vg - (Vv) T Vg. (4) 

Proof. From (3), it follows 

djtp - dtdjip + v k d k djip = dj(d t ip + v k d k ip) - djV k d k <p, 
that is exactly (4). □ 
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3. Generalized Cahn-Hilliard equation 

In this section we introduce the kinetic equation for the phase-field c. As we have remarked in Sect.l, since there is no mass 
transfer from one phase to the other one, the mass of each component is conserved in Q. This means that the evolution of c 
has to be subject to the constraint 

po c(x, t)dv = constant. (5) 



n 

In their original papers, Cahn and Hilliard [7, 8] postulated a generalized mass diffusion equation, valid in the entire two 
phase system, to describe the process of phase separation of two components in a binary alloy under isothermal and isochoric 
conditions. In particular, they assume that the concentration of one of the two metallic components of the alloy obeys the 
equation 

9 t c = V-J, (6) 
where the local diffusion mass flux J satisfies the boundary condition 

J-n| ao = 0. (7) 

In addition, J is assumed to be proportional to the gradient of the generalized chemical potential u,, i.e. 

J = M(c)Vm, 

where M(c) denotes the diffusive mobility and it is a non-negative function eventually depending on the concentration. The 
dependence of mobility on the concentration appears for the first time in the original derivation of the Cahn-Hilliard equation (see 
[7]) and later other authors considered difFerent expressions for M(c) (see for instance [3, 12]). The case M(c) = corresponds 
to a pure transport of the components without diffusion . 
Cahn and Hilliard take p, in the form 

p, = -<yAc+f(c), (8) 

where 7 measures the width of the diffusive layer and f is a double-well function, whose wells represent the two bulk phases. 

Later, this model has been applied to other contexts concerning two phase systems constituted by different substances, for 
instance air and water or oil and water. Furthermore, the Cahn-Hilliard equation has been coupled with a kinetic equation for 
the absolute temperature. For these models thermodynamically consistency and results concerning existence, uniqueness and 
long-time behaviour of the solutions have been proved (see e.g. [5, 18, 26]). 

To account for fluid motion, some authors analyze the so-called Navier-Stokes-Cahn-Hilliard system (see for instance 
[17, 20, 25]). They substitute the partial derivative of c with respect to t with the material derivative in (6), i.e. 

PoC = V • [M(c)V/i] (9) 

where p is given in (8) and they consider a modified Navier-Stokes equation including a surface tension source term for the 
coupling between c and v. In such a way, equation (9) is composed of both a transport term, v • Vc, accounting for mechanical 
effects and due to the presence of the material derivative, and a diffusive term at the right-hand side modelling the chemical 
effects. The drawback of this model is that when slow processes are considered, namely c « d t c, the coupling between c and v 
disappears. 

In our paper we propose a thermodynamically consistent model for phase separation phenomena with a different coupling 
with the fluid motion and including thermal efFects. In particular, we assume that c satisfies equation (9) where the chemical 
potential p is allowed to depend on the curl of the velocity. More precisely, p is taken in the form 

p = -7Ac + e Q F'(c) + [6 + |V x v| 2 ]G'(c), (10) 
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where 7 and 80 are positive constants and F, G are suitable function depending only on c whose expressions are given in the 
sequel. Accordingly, the explicit form of the Cahn-Hillard equation is 

poc = V • {M(c)V [-7AC + 9 F'(c) + (9 + |V x v| 2 ) G'(c)] } . (11) 

We append to such an equation homogeneous Neumann boundary conditions both for the difference concentration and the 
chemical potential, i.e. 

Vc -11130 = 0, VM-n|an = 0. (12) 

The first condition describes a "contact angle" of ir/2 between the diffused interface and the boundary of the domain, while 
the second one means that there is no mass flux through the boundary and it ensures that (5) holds. Indeed, in view of the 
transport and divergence theorems, we have the following equalities: 

d_ 

di 



Pocdv= I pocdv= I M(c)Vp, ■ n da = 0. 

9Q 



A typical choice of the functions F, G is the following: 

F ( c ) = ir-^. G(c) = ^, (13) 



so that when 



is constant, the function 



u = e + |Vxv| 2 (14) 



c 4 c 2 
W(c) = 9 F(c) + uG(c) = 6 — + (u- 6 ) — 



coincides with the double-well function f given in (8). 

Accounting for the explicit expression of F and G, we are able to explain the existence of a critical value for the temperature 
and the curl of the velocity. So, just to this aim, let us neglect the quantity A(— ^Ac) in equation (11) and suppose to fix the 
values of the temperature and the curl of the velocity. Under these hypotheses, the evolution equation for c reads 

p c V • [M(c)W"(c)Vc] = V ■ [K(c)Vc], (15) 

where the diffusivity K is defined as K(c) = M(c)W"(c) and W" is given by 

W"(c) = 39 c 2 + u-e . 



Note that since M(c) is a non-negative function, the qualitative behavior of the solution depends on the function W. Precisely, 
when 9 or |V x v| 2 are sufficiently large (that is u is sufficiently large), the diffusion coefficient K(c) is positive since W is convex 
and it attains a (unique) minimum at c = 0, i.e. the mixed phase is stable. On the other hand, if u < 80, then W has two minima 
at c — ±y/ ^jjp and a local maximum at c = 0. This means that when u < 9o K(c) is negative in the so-called spinodal interval 
(— ci, ci) with 




and it is positive when c < — ci or c > ci. So the mixed phase is unstable and phase separation occurs (see Fig. 1). The constant 
value 9o can be interpreted as the critical temperature of the mixture. 

Notice that equation (15) allows backward and forward diffusion and the corresponding initial problem is classically not well- 
posed from the mathematical point of view. For this reason, in the Cahn-Hilliard equation the additional term A(— "(Ac), which 
accounts for the interfacial energy, appears. 
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W(c) W(c) 




c 



Figure 1. On the left the graphic representation of W(c) when u > Bq. The minimum in zero implies that, when the temperature or the mixing velocity are 
sufficiently large, the mixed phase is stable. 

On the right the graphic representation of W(c) when u < 6q. In this case, we observe two minima and the local maximum corresponding to c = 0: the mixed 
face is unstable. 



4. Governing Equations 

The expression of the chemical potential p. involves the fluid velocity and the temperature. Accordingly, we need to write the 
kinetic equations for these variables. 

By Eq. (1), we are modeling the fluid as an incompressible material. Then the continuity equation provides 

V-v = 0. (16) 

The linear momentum balance equation is taken in the classical form of continuum mechanics, namely 

p v = V ■ T + p b, 

where T is the Cauchy stress tensor and b is the body force (per unit mass). We write T as the sum of three second order 
tensors, i.e. 

T = T" + t + S. 

The first term is related to the classical Cauchy stress tensor in the Navier-Stokes equation, that is 

T = -pi + £/(c)[Vv + (Vv) T ] = -pi + 2u(c)D, 

where p is the pressure (which is a priori unknown since we have supposed the fluid incompressible), 1 is the second order identity 
tensor, u is the viscosity of the fluid and D is the symmetrical part of the gradient of velocity. We stress that u depends on the 
concentration c, so that when c = 1 (or c = —1) v coincides with the viscosity of the fluid A (or B). In the next section, we 
will prove that u(c) > as a consequence of the second law of thermodynamics. 

We add to the usual (symmetric) tensor T" the extra reactive stress T associated with the presence of concentration gradient 
which models the capillary forces due to surface tension (Ericksen's stress, see [13]), i.e. 

t = -7PoVc ® Vc, 

where the parameter 7 is assumed to be positive and it is related to the thickness of the interfacial region. This term occurs 
even on other papers concerning with the Navier-Stokes-Cahn-Hilliard equation (see [19, 25]). 
Finally, we introduce the skew tensor S whose components are defined as 

S,j = £ IJk poG(c)(V x v)k, 

where e^k denotes the Levi-Civita symbol and summation is implied by index repetition. This term makes the tensor T non- 
symmetric, consistent with the presence of internal structure due to the mixture (see [11]). However, S disappears when G = 0, 
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that is in the bulk phases. 

By evaluating the /— th component of the divergence of S, we obtain 

(V • S); = djSij = dj[e uk p G(c)(V x v) k ] = e Uk dj[poG(c)(V x v) k ] 

namely, 

V • S = p V x [G(c)V x v]. 
Accordingly, the velocity v satisfies the modified Navier-Stokes equation 

p v = -Vp + V ■ {2v(c)D} - 7p V • (Vc ® Vc) + p V x [G(c)V X v] + p b. (17) 

To this equation we associate the usual no-slip boundary condition: 

v|an = 0. 

In order to obtain the kinetic equation for the temperature, let us consider the first law of thermodynamics as in [15] or [16] 

p E = V i m + V' c + p h, (18) 

where E is the total energy, V' m ,V' c are respectively the internal mechanical and chemical power whose expressions are given in 
the sequel, and h stands for the rate at which the heat is absorbed by the material. Denoting with T = iv 2 the kinetic energy 
and e the internal energy, which we suppose function of the state a = (6, c, Vc) of the system, we write E = T + e. 

By multiplying equation (17) by v and accounting for (16), we obtain the power balance related to the velocity v, that is 

V = V e 

' m "mi 

with 

V' m = ipo-^v 2 + ^(c)|Vv| 2 + ^(c)(Vv) r : Vv-7p (Vc®Vc) : Vv (19) 
-p G(c)|V x v| 2 , 

Vm = V • [-pv + ^(c)Dv -7p (Vc® Vc)v + p G(c)(V x v) x v] + p b ■ v. (20) 

Similarly, multiplying equation (10) by poc and taking (9) into account, we obtain the power balance related to the 
concentration c, that is 

V' c = VI, 

where 

V' c = p e F(c)+p G{c)[e+\V x v| 2 ] +po7Vc- Vc + M(c)|Vm| 2 , (21) 
VI = V-[po7cVc+M(c)mVm]- (22) 

By means of Lemma 1 that applied to the concentration c yields 

VE ■ Vc = [Vc - (Vv) T Vc] • Vc = Vc ■ Vc - (Vc ® Vc) : Vv, (23) 
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and summing up V' m and V' c , we obtain 



I v 2 + eoF(c)+ I 7 | V c| 2 



+ u{c)\ Vv| 2 + i/(c)(Vv) T : Vv + p o 0G(c) + M(c)\Vp,\ 2 . (24) 



We stress that the fourth term in the rhs of equation (19) is exactly T : Vv, hence by (23), the contribute of T to the internal 
power is enclosed into Po^(^7| Vc| 2 ). 

Moreover, (24) suggests to define the internal energy e as 

e(CT) = eo(0)+0oF(c) + i7|Vc| 2 , (25) 

where eo is a suitable function depending only on the temperature. Then, the total energy E is given by 

E = T + e = iv 2 + e (<9) + 8 F(c) + ^7|Vc| 2 . 

As a consequence, a comparison with (18) yields 

poh = p Q ee' (e) - m(c)|Vv| 2 - i/(c)(Vv) f : Vv - p 9G{c) - M(c)\Vp,\ 2 . (26) 

In our model the Fourier theory of heat conduction will not be modified. Accordingly, the constitutive equation relating the 
heat flux q to the gradient of the temperature, assumes the classical form 

q = -/c(0)V0, (27) 

where k(0) > denotes the thermal conductivity and it depends on the absolute temperature. 
As well known (see e.g. [15]), the thermal balance law is expressed by the following equation 

p h = -V ■ q + p r. 

Comparing with (26), we have 

p ee' (8) - m(c)|Vv| 2 - Kc)(Vv) T : Vv - p o 0G(c) - M(c)|Vm| 2 = -V • q + p r. 
Finally, collecting the equations of motion we write the system of equations: 
V -v =0 

p v = -Vp + V ■ Kc)D] -7p V ■ (Vc® Vc) + p V x [G(c)V X v] +p b 

Poc = V- [M(c)V(-7Ac + e F'(c)+ [e+|Vx v| 2 ] G'(c))] 

p ee' {e) = m(c)|Vv| 2 + m(c)(Vv) t : \7v + p 9G(c) + M(c)\Vp,\ 2 - V ■ q + p r 

and we associate to (28) the boundary conditions 

Vc-n|gn = 0, Vp. • n| S n = 0, 
v|sn = 0, Ve-n| sn = 0. 

and the initial data 

c(x,0) = c (x), v(x, 0) = v (x), 0(x,O) = o (x). 



(28) 
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5. Thermodynamics 

In this section we show that our model is consistent with the second law of thermodynamics written in the Clausius-Duhem 
form. 

Second law of thermodynamics. There exists a function 77, called entropy function, such that 

P0V>-V-(!) + P -f, (29) 



where q is the heat flux vector and r is the external heat supply density. 
We introduce the Helmholtz free energy density ip defined as 

ii(a) = e(ff) - 9t,{g). 

Now we are in a position to prove the following result. 

Theorem 1 The functions q, ip and r\ are compatible with the second law of thermodynamics if and only if the viscosity u, the 
mobility M and the thermal conductivity k are non-negative functions and the free energy ip satisfies the following conditions: 

dgip = - v , d c ti> = e F'(c) + eG'(c), dvctp = iVc. (30) 

Proof. In order to obtain compatibility with thermodynamics, we have to prove that system (28) with constitutive equations 
(26), (27) agrees with inequality (29), which, by means of the thermal balance law, becomes 

PoT?e > ^Ve -q + po/7. (31) 

Moreover, by the free energy density ip = e — 9r), (31) can be written as 

potp - poe + p 6ri + p h + \v9 • q < 0. 

In view of (25)-(26) we have 

po {d e ip + r?) 9 + po [dcip - 9 F'(c) - 9G'(c)] c + p {d Vc ip - tVc) • Vc 

-u (c)|Vv| 2 - t^(c)(Vv) T : Vv - M(c)|Vm| 2 + ^V(9 ■ q < 0. (32) 
From definition (10), it follows that 

Vp, = -7V(Ac) + V{e F'(c) + [9 + |V x v| 2 ]G'(c)}. 

Since V(Ac) may be chosen arbitrarily, V/j. may be chosen arbitrarily too. Accordingly, by standard arguments, we deduce the 
following conditions: 

datl> = -v, d c ^ = 9 Q F'{c) + 9G'{c), cVcV = 7V c 
and, in view of (27), we write inequality (32) in the form 

-i^(c)|Vv| 2 - Kc)(Vv) T : Vv - M(c)|Vm] 2 - 4r^|V0| 2 < 0. 



In addition, relation 



(Vv) T : Vv = tr(|Vv| 2 ) > 0, 
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allows us to conclude that the function v(c), M(c) and k{9) are non-negative. □ 
From (30) it follows that the free energy density tp and the entropy 77 are given, up to a constant, as 

ip = 6 F(c) + 6G(c) + 'l\X7c\ 2 + Md), (33) 
r? = _ 9 ^ = _G(c)-^(0), (34) 

where ipo is a suitable function (depending only on 9) which ensures the validity of the condition ip = e — 776 = e + dgip9. A 
substitution of (25) and (33) leads to the equality 

Me) = e (e) + ij' (e)e. 

Thus, ipo is given by 
with C > and 

In particular, if we let eo = C6, where C denotes the specific heat, we recover the standard form of ipo and 77, i.e. 

ipo =C9(1 -\n9), n = -G(c)-C\nO. 
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